Modelling the structure of clusters of Ceo molecules 
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We locate putative global minima for (C6o)jv clusters modelled by the potential of Pacheco and 
Prates-Ramalho up to 7V=105. These minima are based on icosahedral packing up to A'^=15, but 
above this size the lowest-energy structures are decahedral or close-packed. Although structures 
based on the 98-molecule Leary tetrahedron, which have been inferred from experiment, are 
not lowest in energy for this potential, an examination of the energetics of a growth sequence 
leading to the Leary tetrahedron lends further support to the experimental assignments. An 
analysis of the potential energy surface topography and the thermodynamics of two example 
clusters indicates that the multiple-funnel topography is likely to have a strong influence on 
the dynamics of structure formation and that solid-solid transitions driven by differences in 
vibrational entropy are likely to be common. 
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I. INTRODUCTION 

There has been much interest in the condensed- 
phase properties of Cgo molecules because of their un- 
usual intermolecular potential. At high temperature 
the molecules can rotate freely and so they act as 
large 'pseudo-atoms' with interactions that are extremely 
short-ranged relative to the large equilibrium pair sepa- 
ration. Consequently the properties of Ceo can extend 
beyond the range of behaviour observed for atomic mate- 
rials. For example, there aralheoretical prediction|.tl. 
the liquid phase is unstableEla or marginally stabledui 
with the precise results depencling on the potential (and 
somewhat on the methodologyCl) used. 

There has also been considerable interest in the struc- 
tural properties of clusters of Ceo molecules, prompted 
by the first experiments of Martin et al. on positively- 
charged clusters. H This mass spectroscopic study revealed 
magic numbers for clusters with less than 150 molecules 
that are injiicative of structures based upon Mackay 
icosahedraE The stability of the icosahedral (dso^p, 
has since been further illustrated by Hansen et a/..ll3'til 
However, subsequent calculations using the spherically- 
averaged Girifalco potential^ found that icosah£dX|al 
structures are only lowest in energy up to A^=13JljO'E£l 
Above this size the structures are either decahedral or 
close-packed. The rapid emergence of bulk-like struc- 
tures for this potential contrasts with many other atomic 



clusters where icosahedral structures can persist up to 
large sizes, for example, up to at least 20 000 atoms for 
sodium.llj This behaviour, and the marginal stability of 
the liquid phase, has a common origin in the narrow- 
ness of the Ceo intermolecular potential well. As the 
range of a potential is decreased, there is an increas- 
ing energetic penalty for strained stcuctures, be they 
icosahedraHl or liquid configurations Jl30 which involve 
nearest-neighbour distances that deviate from the equi- 
librium pair distance. 

One possible cause of the discrepancy between experi- 
ment and theoretical calculations is the isotropic nature 
of the Girifalco potential. However, calculations using an 
all-atom model also found that icosahedra are only low- 
est in energy at small sizes, albeit to sligktJKpLarger sizes 
(iV=16) than for the Girifalco potentialMBB 

This situation led Shvartsburg and Jarrold to inves- 
tigate the possibility of distinguishing icosahedral from 
decahedral and close-packed isomers by their mobil- 
ity. However, the differences in the computed mobili- 
ties are sniaJl and were comparable to the experimental 
resolution.L^ More recently, Branz et al. performed some 
new mass spectroscopic experiments on (Ceo)^ clusters 
to try to obtain further insights into the difference be- 
tween theory and experiment .Ej They found that the ob- 
served structures are independent of the sign and mag- 
nitude of the charge. However, temperature was found 
to be a key variable in determining the structure. The 
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initial cold as-grown clusters showed no magic numbers. 
Only on annealing at higher temperatures are the magic 
numbers revealed, as the relative evaporation rates cause 
larger populations to develop in those clusters that are 
more resistant to evaporation. After annealing at 490 K, 
as in the previous experiment, magic numbers consis- 
tent with icosahedral clusters are obtained. However, 
annealing at 585 K reveals a new set of magic Jium- 
bers that correspond to non- icosahedral clusters Bj As 
well as sizes associated with face-centred-cubic (fee) and 
decahedral packing, particularly prominent were magic 
numbers associated with the recently discovered Leary 
tetrahedron.cj 

There are two possible explanations for this tempera- 
ture dependence. Firstly, the results could reflect changes 
in the thermodynamically most stable structures with 
temperature. Such solid-solid t r aiia itio i lis . . ha^^e | now been 
observed in a variety of systems. EZlEZrESnjI^d'EllLSo How- 
ever, for all these examples icosahedral structures are 
more stable at higher temperature, and are faicpured be- 
cause they have a larger vibrational entropy.o Further- 
more, the theoretical calculations for (C6o)iv suggest that 
the non-icosahedral structures are lowest in energy and 
so would be favoured at low temperature. Indeed, such 
a transition to a high-temperature icosahedral structure 
has been located for (C6o)i4 interacting with the Giri- 
falco potentialH. and has also been suggested for some 
larger clusters. O As the behaviour of these solid-solid 
transitions is opposite to that of the experiment, this ex- 
planation is unlikely. 

Secondly, the icosahedral structures could be more ac- 
cessible during the growth and annealing of the clusters, 
with escape from the icosahedral region of configuration 
space into the state with lowest free energy only being 
possible at the highest temperatures. This interpreta- 
tion has support from a variety of sources. For exam- 
ple, it has been shown that it is possible, under certain 
conditions, to preferentially grow metastable icosahedral 
clusters for silver Furthermore, analyses of the poten- 
tial energy surface (PES) topography of those Lennard- 
Jones clusters with non-icosahedral global minima has 
shown that large energy barriers exist for interconversion 
of the icosahedral isomers and the global |»inimum, and 
that the icosahedral funnel is much wider .E3 Thus, on re- 
laxation down the PES, the system is likely to become 
trapped in the icosahedral region of configuration space. 
This trapping is partly due to the greater structural simi- 
larity to the polytetrahedral configurations typical of the 
liquid, but also to the general energetic preference that 
Lennard-Jones clusters have for icosahedral geometries. Ell 
The latter would not hold for clusters of Ceo molecules. 
Finally, short-ranged potentials have been shown to lead 
to relatively large_barriersE§E3 and consequently to much 
slower dynamicso Therefore, it is plausible that kinetic 
effects could dominate for (C6o)Ar clusters on the exper- 
imental time scales, even for the relatively small clusters 
studied. 

Many of the studies of the condensed-phase properties 
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FIG. 1: Comparison of the PPR and Girifalco pair po- 
tentials. 

of Ceo have used the Girifalco potential.liliiillillilli.El 

However, a new single-site potential has recently been 
derived by Pacheco and Prates-Ramalho (PPR) that 
gives an improved description of many properties of bulk 
Ceo, including, for example, thii. dependence of the den- 
sity of the crystal on pressures Furthermore, the high- 
temperature thermodynamics of clusters interacting with 
this potential are in ijeasonable agreement with those of 
an all-atom potential.c2l Here, we see if this potential can 
hel|P,to explain the experimental observations of Branz et 
aZ.Bj Global optimizatimiias been previously attempted 
for the PPR potentiaLEJeJ However, the conclusions of 
these papers are contradictory and, as we shall see, many 
of the putative global minima are sub-optimal. Here, as 
well as locating putative global minima for this poten- 
tial, we also consider the most stable growth sequence for 
structures based on the experimentally observed Leary 
tetrahedra, and examine the PES topography and the 
thermodynamics of two example clusters. 

II. METHODS 

The PPR potential consists of a pair potential plus_a. 
three-body term of the standard Axilrod- Teller form.c£l 
The potential energy of the cluster is therefore given by 

= y^axTirij) + Cat ^ VAT(7'y,fj/c,?'jfe), (1) 

i<j i<j<k 

where Cat gives the magnitude of the Axilrod- Teller 
term. The pair potential consists of a Morse form for the 
short-range repulsion, a van der Waals expansion for the 
long-range attraction and a Fermi function to describe 
the crossover between these two regimes. The forms of 
these functions and the associated parameters are given 
in Ref. H 
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As one can see from Figure |l|, the PPR pair potential 
has a softer repulsion than the Girifalco potential, leading 
to a wider potential well. This effect can be quantified 
by matching the second derivative at the bottom of the 
well to that of the Morse potential, 

VAiir) = eexp [p{l - r/r^^)] (exp [p{l - r/r^^)] - 2) , 

(2) 

where e is the pair well depth and r^q is the equilibrium 
pair distance. The Morse potential becomes increasingly 
narrow as the parameter p increases. The curvature at 
the bottom of the Morse well is 2p^ when the units of 
energy and distance are the pair well depth and equi- 
librium pair distance. This result can then be used to 
obtain a value of pcs, a measure of the effective range, 
for each potential. This analysis has been dotie previ- 
ously for the Girifalco potential: /9^=13.62.c£l For the 
PPR potential, however, p^g^=11.28. This difference in 
Pefi has a wcll-updetstood effect on the resulting proper- 
ties of the liquidllaO'cZl and of clusters .EZlc3 For example, 
it will make the bulk PPR liquid phase more stable, as 
has been observed.H More importantly for this study, it 
will make icosahedral structures more stable than for the 
Girifalco potential. 

The three-body Axilrod- Teller term always gives rise 
to a positive contribution to the energy for a compact 
structure. For example, its inclusion leads ta a 6% in- 
crease in the energy of the bulk Cgo crystal.cZl However, 
in the global optimization study of Ref. the supposed 
global minima for the full PPR potential had a lower po- 
tential energy than those when only the pair potential 
was used. This is clearly wrong and presumably must 
have been due to an error when coding the Axilrod- Teller 
term. Furthermore, in the other optimization study of 
PPR clusters, putative global minima were only located 
for the pair potential, because of rtiie greater computa- 
tional cost of the three-body term.E3 

The most unfavourable common configuration 
for the three-body Axilrod- Teller term is three 
nearest-neighbour molecules arranged in an equilateral 
trianglcEj Therefore, one would expect the three-body 
energy to be larger for those structures that have more 
nearest neighbours, more polytetrahedral character and 
close-packed surfaces. Thus the AxilrodiTeller term 
disfavours icosahedral structures the most.E3 

To locate the global minima for the PPR po- 
tential we used basin-hoppingEll (Monte Carlo plus 
minimizationt3) , which has proved to be a very effec- 
tive mctked- for a variety of cluster sys 
approachEilE^ and the reasons for its succes£ 
been described in detail elsewhere. The only specific 
modification for the present application was to reduce the 
computational cost of the minimizations by only switch- 
ing on the three-body interactions close to convergence. 
Such a 'guiding function' approach iias previously been 
suggested and exploited by Hartke.EJ As well as basin- 
hopping, we also reoptimized a large database of minima 
that we had obtained from previous optimization studies 
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on Lennard-Jonesj^ MorseSS and GirifalcoEl clusters. 
Most of the global minima that we located were contained 
within this database. 

To generate the samples of minima from which the dis- 
connectivity graphs and the thermodynamics in Section 
[V were calculated, we used thft-aame, methods as-|those 
we have previously applied to LJl2Il3EZI and MorseEj clus- 
ters. We thereby obtain large samples of connected min- 
ima and transition states that provide good representa- 
tions of the low-energy regions of the PES. The approach, 
involves repeated applications of eigenvector-following^ 
to find new transition states and theptninima they con- 
nect, as described in detail elsewhere.Ej 



III. GLOBAL MINIMA 

Putative global minima were located for the full PPR 
potential up to A^=105, the size range of interest for com- 
parison, with the high temperature experiments of Branz 
et a/..Cj The energies and point groups of these structures 
are given in Table || and coordinates are available on tka 
world-wide web from the Cambridge Cluster Database.E2l 

As expected, the energies of the global minima dif- 
fer from those reported in Ref. because of the er- 
ror in the Axilrod- Teller term. To facilitate comparisons 
with the other results in Refs. |^ and we also re- 
port the energies of these global minima when reopti- 
mized for the PPR pair potential (Table |). The re- 
sulting energies will not necessarily be the global optima 
for the pair potential, but they should provide a good 
upper bound. Comparisons with these energies show 
that the putative global minima for the PPR pair po- 
tential given in Ref. ^ are sub-optimal for > 19 and 
those in Ref. [34| are sub-optimal for N=18, 26, 33-36 and 
N > 40. Only for N=16 and 29L-a|p, lower two-body en- 
ergies obtained in these papers,OCj indicating that for 
these cluster sizes the introduction of the Axilrod- Teller 
term leads to a change in the structure of the global 
minimum. For A=16 the relevant minimum is icosa- 
hedral with £;=-13.006550eV and £;2=-13.447202 eV, 
and for A=29 it is an alternative decahedral minimum 
with £;=-27.540094eV and £;2=-28. 516282 eV. 

In Figure ^ja we plot the energies of the PPR global 
minima so that particularly stable clusters stand out. 
We also provide a comparison with results for the Giri- 
falco potential (Figure ||b), partly because a significant 
number of the putative global minima in Ref. |l^ have 
now been improved (an up-to-date, list is maintained at 
the Cambridge Cluster DatabaseEJ). A selection of the 
global minima are depicted in Figure ^. 

At iV=13 the icosahedral global minimum is notice- 
ably more stable than for the Girifalco potential, which 
is consistent with the PPR potential's slightly wider well. 
Furthermore, icosahedral structures are lowest in energy 
up to 15 molecules, or up to 16 molecules if the Axilrod- 
Teller term in the potential is not included — as noted in 
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FIG. 3: A selection of the global minima for the PPR potential. Each molecule is represented by a point at the 
molecular centre. The two views of the 48- and 58-molecule clusters illustrate the relationship of these structures to 
both decahedra and the Leary tetrahedron. 
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FIG. 2: Energies of the global minima for (C6o)jv clusters 
modelled by the (a) PPR and (b) Girifalco potentials. 
The energy zero is taken to be i?avo, a four-parameter fit 
to the energies of the global minima. For the PPR poten- 
tial E^^JeY= -1.6537iV + 2.1901iV2/3 + 0.1263iVi/3 - 
0.7467. For the Girifalco potential E^^a/eV= -1.789iV-F 
2.294A^2/3 _,_ o.5907iVi/3 - 1.292. 



Section || this term slightly disfavours icosahedral struc- 
tures. However, contrary to the claims of Refs. ^ and ^ 
there are no further icosahedral global minimum above 
this size. Instead decahedral or close-packed clusters are 
lowest in energy. 

To illustrate how the relative stabilities of icosahedral 
structures with respect to the global minimum evolve 
with size, we give the energy differences for examples 
where particularly stable fee and icosahedral forms are 
available. At A^=38 the difference in energy between the 
fee truncated octahedron and the lowest-energy icosahe- 
dral structure is 0.582 eV. However, at N=55 the Mackay 
icosahedron is 0.297 eV above the decahedral global min- 
imum. When the Axilrod- Teller interactions are not in- 
cluded this difference decreases to 0.159 eV, again illus- 
trating that this term slightly disfavours icosahedra. In 
contrast, the energy difference for the Girifalco potential 
is a considerably larger 1.997 eV. These results are con- 
sistent with the behaviour expected from the respective 
values of p^s: the Mackay icosahedron is lowest in energy 
for the Morse potential up to p=11.15 and then becomes 
increasingly unstable as the width of the potential de- 
creases further. 

The pronounced minima in Figure ^ correspond to par- 
ticularly low-energy structures that could potentially give 
rise to magic numbers. Comparing (a) and (b) of Figure |^ 
it is apparent that the patterns are very similar, indicat- 



FIG. 4: Particularly stable clusters on the growth sequence leading to the Leary tetrahedron. When two viewpoints 
are given they correspond to the front and back of the cluster. 



ing that the two Ceo potentials favour similar structures. 
Prominent in both figures are the same particularly sta- 
ble decahedral and close-packed forms, most of which are 
illustrated in Figure ^. The main differences are that, 
for the Girifalco potential, the feature due to (Ceo) 13 is 
much less pronounced, for the reasons explained above, 
and that the close-packed structures are somewhat more 
favoured. The latter feature is again due to the shorter 
range of the Girifalco potential, as is the greater number 
of close-packed global minima for the Girifalco potential. 

There is some correspondence between Figure 2 
and the high —temperature magic numbers observed 
experimentally£j Small peaks in the mass spectrum at 
N—31, 33 and 35 probably correspond to a series of 
small decahedra. Similarly, the peak at iV=38 can be 
assigned to the fee truncated octahedron, and the peaks 
at N=6i, 71 and 75 are probably due to Marks deca- 
hedra. Particularly interesting is the feature in Figure |^ 
corresponding to the decahedral (Ceo)48 global minimum, 
because the most prominent peak in the mass spectrum 
occurs at A^=48. Like a number of the other decahe- 
dral global minima, this structure has an overlayer on 
the {111} faces in sites which are hep with respect to the 
five slightly strained fee tetrahedra that are the basis of 
decahedral structures. In fact, many of these structures 
are fragments of larger Mackay icosahedra with the apex 
of the decahedron corresponding to what would become 
the centre of the icosahedron. Indeed, growth simula- 
tions for silver have indicated that the addition of this 
overlayer does provide a natural pathway to the growth 
of larger icosahedra.E3 

However, many of the particularly stable structures, 
especially for A'^>50, do not correspond to magic numbers 
in the experiment. There are no experimental peaks cor- 
responding to the close-packed tetrahedra at A^=59 and 
100, nor to the twinned truncated octahedra at A^=50 
and 79 and the 101-molecule Marks decahedron. Instead, 
peaks occur in the mass spectrum at iV=58, 61, 68, 77, 
84, 91, 96 and 98, which cannot be simply explained by 



reference to the global minima of model Cen^potentials. 
These peaks have previously been assignedcfl tq-atruc- 
tures based on a 98-molecule Leary tetrahedron.E£l This 
recently discovered structure is illustrated in Figure ^. It 
can be described in terms of five fee tetrahedra arranged 
into a stellated tetrahedron with apex molecules removed 
and 7-molecule hexagonal patches covering the edges of 
the central tetrahedron. That this structure is not the 
global minimum for the PPR or Girifalco potential is con- 
sistent with the values of PeS- For the Morse potential it 
is only the global minimum for 6.91 < p < 9.45 because 
it has a strain energy intermediate between decahedra 
and icosahedra. 

To attempt to add further weight to the experimen- 
tal assignments we calculated the energies of a series of 
structures derived from the Leary tetrahedron. Although 
they are not lowest in energy, the size variation of their 
energies, plotted for 60 < A'' < 100 in Figure ||, is insight- 
ful. There is remarkable agreement between the features 
in this graph and the remaining magic numbers. 

The stable structures for the Leary tetrahedral growth 
sequence are illustrated in Figure |. The 91-molecule 
structure is derived by removing a hexagonal patch from 
over one of the edges of the central tetrahedron. In 
contrast, to obtain the 89- molecule structure one of the 
points of the Leary tetrahedron is further truncated. By 
a similar truncation the {C6o)9i global minimum can be 
derived from the 100-molecule tetrahedral global mini- 
mum (Figure ^). Most of the other structures can be 
obtained by a combination of these two changes. The 
84-molecule structure is derived by a truncation and the 
removal of an adjacent patch. Similarly, truncation of 
two of the fee tetrahedra and the removal of the common 
patch gives the 77-molecule structure. The 68-molecule 
structure is derived by completely removing one of the fee 
tetrahedra and adjacent patches, thus giving a structure 
that is also a fragment of the 147-molecule Mackay icosa- 
hedron. The truncation of a further tetrahedron leads to 
the 61-molecule structure. 



6 




60 65 70 75 80 85 90 95 



N 

FIG. 5: Energies of structures based upon the 98- 
molecule Leary tetrahedron (sohd hue) compared to the 
energies of the global minima (dashed line) for the PPR 
potential. 



Interestingly, the latter two structures are global 
minima for Morse clusters, but were not-, located 
in the most recent optimization study.cS They 
have energies of £;6i(/o=10)=-252.488332 e and 
-£'68(p=10)=— 286.643320 e, and are lowest in energy in 
the ranges 9.42 < p < 10.34 and 7.40 < p < 11.55, 
respectively. 

The prominence of the peak in the mass spectrum at 
7V=48 also fits with this preference for structures based 
upon the Leary tetrahedron. The second view of the dec- 
ahedral (C6o)48 global minimum in Figure ^ shows that it 
is also a fragment of the Leary tetrahedron. By adding an 
identical overlayer to the bottom {111} faces of this clus- 
ter the decahedral (C6o)58 global minimum is obtained, 
which is again a fragment of the Leary tetrahedron. Al- 
though this cluster is not especially stable compared to 
nearby sizes for the PPR potential, its relationship to 
the Leary tetrahedron provides a basis for confidently 
assigning the remaining magic number at N=58 to this 
structure. 



IV. (C6o)38 AND (C6o)55 

Although the location of the global minima should be 
the first step in the theoretical determination of a clufe 
ter's structure, the thermodynamics^ and dynamicsEll 
can also play an important role. Here, we examine the 
behaviour of two example clusters in more detail. We 
choose (C6o)38 and (€60)55 because particularly stable 
fee (38) and icosahedral (55) structures are possible. Fur- 
thermore, these sizes have been extensivebi,jS±juiiied_J'pi: 
the longer-ranged Lennard- Jones potential.L3ElE3c3l£3S 

In order to construct disconnectivity graphs for these 
two clusters we have generated large samples of min- 



ima and transition states, plus the pathways connect- 
ing them. As this is a computationally demanding task 
only the two-body component of the PPR potential was 
used. The graphs for the full potential would be very 
similar; the main effect of the introduction of the three- 
body term would be to displace the icosahedral regions 
of configuration space further up in energy. We located 

38 558 minima and 39 959 transition states for A^=38 and 

39 043 minima and 39 845_t|i:ansition states for A^=55. 
Disconnectivity graphsoLj provide apcepijesentation of 

the barriers between minima on a PES.OEJ In a discon- 
nectivity graph, each line ends at the energy of a mini- 
mum. At a series of equally-spaced energy levels we de- 
termine which (sets of) minima are connected by paths 
that never exceed that energy. We then join up the lines 
in the disconnectivity graph at the energy level where the 
corresponding (sets of) minima first become c onpg cted. 
In a disconnectivity graph an ideal single- funneQEZI PES 
would be represented by a single dominant stem associ- 
ated with the global minimum to which the other min- 
ima directly join. For a multiple-funnel PES there would 
be a number of major stems that only join at high en- 
ergy. A single funnel PES is typically associated with 
efficient relaxation to the global minimum, whereas for a 
multiple-funnel PES there is a separation of time scales 
between relaxation-t|0 a low-energy structure and inter- 
funnel equilibrium 

For (C6o)38 and (€50)55 it is not possible to charac- 
terize the whole PES because of the huge numbers of 
minima. Besides, even if we could perform such a char- 
acterization, any attempt at representation using a dis- 
connectivity graph would just be obscured by the density 
of lines. Therefore, we concentrate on the low-energy re- 
gions of the PES that are of most importance when con- 
sidering structure, and only include the lines leading to 
the lowest 250 minima. 

For both clusters the PES tography has a multiple- 
funnel character, where each funnel corresponds to a dif- 
ferent structural type (Figure |^) . For this short-ranged 
potential, there are a number of competing low-energy 
morphologies that are separated by large energetic bar- 
riers. As the barriers between minima of the same struc- 
tural type are typically much smaller, the graphs clearly 
separate the different morphologies. 

This behaviour contrasts with that for Lennard- Jones 
clusters, which (with some notable exceptions) typically 
have a single-funnel topography that is dominated by 
icosahedral structures£3 However, for the PPR potential 
the difference in energy between close-packed and deca- 
hedral structures is small. Furthermore, there are often 
a number of close-packed forms possible that have sig- 
nificant structural differences. For example, for (€50)55 
the close-packed region of configuration space divides up 
into structures that are based on the 50-molecule global 
minimum (the twinned truncated octahedron) and those 
based on the 59-molecule tetrahedral global minimum 
(Figure ^). Additionally, as well as a funnel leading to 
the conventional decahedral global minimum, there is a 
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FIG. 6: Disconnectivity graphs for (a) (C6o)38 and (b) (C6o)55- Lines corresponding to the lowest 250 minima are 
represented. Structural labels have been placed adjacent to the lines corresponding to the different funnels on the 
PES, where, as in Table ||, I stands for icosahedral, D for decahedral (D+ have an overlayer on the {111} faces), F 
for fee and C for close-packed. The numbers refer to the number of minima below the corresponding node in the 
full graph. Pictures of the lowest-energy minima for the different structural types have been placed adjacent to the 
corresponding lines. 



low-energy region of configuration space corresponding 
to decahedral structures that have the {111} faces par- 
tially covered. These latter structures can be constructed 
by adding molecules to the 48-molecule global minimum. 
This complexity of the (Ceo) 55 disconnectivity graph con- 
trasts with that of (Ceo) 38: for which there are only fun- 
nels associated with the three basic morphologies. 

The disconnectivity graphs can tell us a lot about the 
dynamics of these two clusters. On the disconnectivity 
graphs we give the numbers of minima in our sample 
that are associated with each funnel. Although, these 
are large underestimates of the true values because of 
the incompleteness of our samples, they do provide a rea- 
sonable indication of the relative numbers of low-energy 
minima of each structural type, and hence of the widths 
of the funnels in configurational terms. It is particularly 



noticeable that there are far fewer low-energy icosahe- 
dral minima, as expected for a short-ranged potential. 
By contrast, the number of decahedral and close-packed 
structures is of the same order. 

Based on this information, on relaxation down the PES 
one would expect the system to be more likely to end up 
in a close-packed or decahedral configuration. However, 
there are other factors that can influence the dynamical 
behaviour. There is a vibrational, as well as a config- 
urational, contribution to the width of a funnel. This 
term is larger for the icosahedral funnel because the as- 
sociated minima are generally flatter and wider. Fur- 
thermore, it may also be that the icosahedral structures 
lie, in some sense, closer to the liquid-like region of con- 
figuration. space because of their greater polytetrahedral 
character.E3C3 
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Structure formation in clusters does not necessarily 
occur by relaxation from a high-energy state, but can 
also occur by growth from a smaller pre-existing solid 
clustcr.EJ Although the width of the funnels is less im- 
portant in this case, the disconncctivity graphs can still 
be useful because they give an idea of the kinetic stabil- 
ity of the possible structures. The large barriers between 
funnels suggest that, except at high temperature, trap- 
ping will occur within the funnel of the existing structure, 
even if the structure is metastable. Indeed escape from 
a funnel will become more difficult with increasing size 
because the intcrfunncl barriers become larger. 

Therefore, during cluster growth, except at high tem- 
peratures and extremely long time scales, the smaller 
clusters serves as templates for growth and only the op- 
timization of the position of additional loosely-bound 
molecules is likely to occur. If icosahedral structures are 
preferred at the last size for which equilibrium is pos- 
sible, then larger icosahedra will result. Thus, the dis- 
conncctivity graphs give additional support to the idea 
that icosahedra are observed experimentally because of 
kinetic trapping at low temperature. 

To examine the low-temperature thermodynamics of 
the two clusters Hirqpestion we use the harmonic super- 
position method.EaEj In this approach the partition func- 
tion is constructed by summing the partition functions of 
the individual minima on the PES. In the harmonic ap- 
proximation, we obtain 

^ ^ n,e^p{-f3EjkT) 

where Ei is the energy of minimum i, Vi is the geomet- 
ric mean frequency and rii is the number of permuta- 
tional isomers {2Nl/hi, where hi is the order of the point 
group). 

This approach is unable to reproduce high tempera- 
ture behaviour such as melting in the present cases, be- 
cause no effort has been made to compensate for the in- 
completeness of our samples of minima, thus leading to 
an underestimate of the configurational entropy of the 
high-energy states. It is possible to overcome this lim- 
itation by n si Tig information obtained from an ergodic 
simulation,HEj but this is not attempted because ergod- 
icity is hard to achieve for these clusters and because we 
are more interested in solid-solid transitions, for which 
we only need to characterize the low-energy regions of 
the PES. Instead we used parallel-tempering to locate 
the point at which the cluster melts or boils. It iSrsJsei 
possible to introduce anharmonicity into this scheme,E3'Lil 
however this additional level of complexity is unnecessary 
for the insights we are seeking to achieve here. 

The advantage of the superposition method is the ease 
with which the low-temperature thermodynamics can be 
obtained even when, as with (Ceo)^ clusters, there are 
large interfunnel energy barriers present. In fact it is not 
clear if any other method would be able to probe this 
regime, because of the difficulty of achieving ergodicity. 
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FIG. 7: Heat capacity curves for (a) (C6o)38 and (050)55 
obtained from our samples of minima using the harmonic 
superposition method. The solid lines are for the full po- 
tential and the dashed lines are for the two-body poten- 
tial. 



Even techniques based on parallel tempering,[l30 which 
are the only direct simulation methods that have been 
shown to succeed for challenging cases invahiing low- 
temperature solid-solid transitions in clusters ,E/Lj failed 
for (060)55. The success of these methods relies on the 
coupling of the low-temperature runs to ergodic high- 
temperature runs involving molten clusters. However, 
because of the narrowness of any range of stability for 
the liquid-like state of clusters of Oeo moleculescirCJ this 
is very hard to achieve. 

Heat capacity curves for the two clusters are presented 
in Figure |7[ These were calculated using the same sam- 
ples of minima as for the disconncctivity graphs. It was 
also possible to reoptimize all the minima for the full po- 
tential. Both clusters show some evidence of solid-solid 
transitions. 

For (000)55 the energy gap between the icosahedral 
minima and the global minimum is relatively small and 
so a transition to the Mackay icosahedron occurs sig- 
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nificantly below the melting temperature. For the full 
potential this transition occurs at 377 K and, consistent 
with the reduction in the energy gap, at 184 K when the 
three-body term is neglected. Interestingly, there is also 
a small peak at '^45 K, which corresponds to a redistri- 
bution of the equilibrium occupation probability amongst 
the five lowest-energy decahedral isomers. These minima 
are almost isoenergetic and correspond to the five differ- 
ent positions that the four-coordinate surface molecules 
can occupy. The high-temperature peak corresponds to 
the beginnings of the melting/boiling transition, which 
as expected is underestimated by the harmonic super- 
position method. Parallel tempering indicates that the 
actual peak for this transition is much larger and occurs 
at ~950K, although the result is somewhat dependent 
on the size of the constraining box in which the cluster 
is placed. 

For (C6o)38 the calculated heat capacity curve is much 
simpler and only has a single peak. This corresponds 
to the transition out of the truncated octahedron, first 
roughly equally into the low-energy decahedral minima 
and lowest-energy icosahedral minimum, and then at 
slightly higher temperature into the higher-energy states. 
As these transitions are so close together in temperature 
there is only a single peak in the heat capacity, which is 
centred on the initial transition because of the underes- 
timation of the latent heat of melting by the present ap- 
proach. Parallel tempering simulations suggest that the 
preliminary solid-solid transition gives rise to a shoulder 
in the melting peak, which occurs at ~900K. 

Although we have only considered two sizes in de- 
tail, the large temperature window for which the Mackay 
icosahedron is most stable indicates that solid-solid tran- 
sitions are likely to be common for other sizes. The 
results confirm what has been emphasized elsewhere, 
namely that crossover sizes at which theudominant struc- 
ture changes depend on temperature.Lj So, although 
none of the clusters above iV=15 have icosahedral global 
minima, icosahedra can still be thermodynamically most 
stable for some range of temperature above this size. For 
example, similar calculations indicate that for (C6o)i9 
icosahedral structures are favoured above 333 K. How- 
ever, these results cannot explain the experimental ob- 
servation of icosahedra, because they are seen at low not 
high temperature. 

As for the solid-solid transitions that hayejbeen investi- 
gated in detail for Lennard- Jones clusters ,£21113 the above 
transitions are driven by the larger vibrational entropy 
of the icosahedral structures. For A^=38 'S'icos:i7dcca:i^cp= 
0.920:0.977:1 and for iV=:55 0.909:1:0.999, where we 
have used the values for the lowest-energy minimum of 
each structural type. These differences are much larger 
than for Lennard- Jones clusters, where, for example, 
i'icos:i^dcca:i^cp=l:l-000:0.990 for LJ55 and the systematic 
differences between compact sequences of structures are 
no larger thaa 2% (Mackay icosahedra having the lowest 
frequencies) Jl3 

F is a measure of how the energy increases as the struc- 
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FIG. 8: The dependence of v on i?strain for (C6o)38 
minima with nnn=139 (crosses) and (C6o)55 minima with 
?inn=217 (diamonds). These samples contain 8408 and 
7337 minima, respectively, and correspond to the largest 
subsets of the complete samples that have the same value 
of rinn. The mean frequency is measured with respect to 
that for the global minimum. 



ture is distorted, and so we can gain insight into the ori- 
gins of the dependence of v on structure if we examine the 
different contributions to the energy. For a pair potential 

E = — nnnfi + -^strain + £^nnn, (4) 

where rinn is the number of nearest neighbours (defined 
using a distance criterion ?'o), the strain energy, the ener- 
getic penalty for distances deviating from the equilibrium 
pair distance, is given by 

i<j,rij<ro 

and the energy of non-nearest neighbours by 

Ennn^ ^ V{n,). (6) 
i<j,rij>ro 

As Ennn IS Small for a short-ranged potential and does 
not vary strongly with distance the main contributions 
to V come from nearest neighbours. Therefore, it is clear 
that the more nearest neighbours a structure has the 
greater the effect on the energy. Thus, compact struc- 
tures with few low-coordinate surface molecules have 
larger vibrational frequencies. This is the reason that V 
generally decreases as the energy of the minima increases. 

The strain energy also has a significant effect on the 
vibrational frequencies. For an unstrained structure all 
nearest neighbours lie at the equilibrium separation, and 
so any distortion of the structure leads to an increase in 
all the individual nearest-neighbour pair energies. For 
a strained structure, there is a dispersion of nearest- 
neighbour distances and so a distortion is likely to cause 
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some of them to deviate further from the equihbrium dis- 
tance, but others to become closer. This compensation 
leads to a smaller rise in the energy, and hence to smaller 
vibrational frequencies. When the effect is isolated by 
only comparing minima with the same number of near- 
est neighbours, this trend is clear (Figure 

As the curvature of the pair potential increases, the 
effect of strain becomes larger. We can therefore explain 
both the much lower frequencies for the icosahedral struc- 
tures in our two examples and why this effect is more dra- 
matic than for Lennard-Jones clusters. Therefore, as the 
range of the potential is decreased, the increasing energy 
gap between icosahedra and the lowest-energy structures 
is compensated by a lower relative mean vibrational fre- 
quency and so it is still possible that icosahedral struc- 
tures are stable for some range of temperature, even when 
the potential is relatively short-ranged. 

The differences in V between close-packed and decahe- 
dral structures are smaller, because of the smaller dif- 
ferences in the strain energy. Indeed, for {C6o)55 the 
lowest-energy decahedral minima actually has a higher V 
because the effect of a larger i^strain is more than com- 
pensated by the effect of a greater n,in- 

As Leary tetrahedra have a strain energy intermedi- 
ate between icosahedra and decahedra, one would ex- 
pect that their mean vibrational frequency is lower than 
for decahedral and close-packed clusters. Indeed, this is 
the case. For example, for (C6o)98 Vicos^i.c&ry'^dcc&^cp= 
0.925:0.972:0.997:1. Although the Leary tetrahedron 
does develop a small equilibrium population near to the 
melting point in the current model, there is no clear tran- 
sition at which the Leary tetrahedron becomes the dom- 
inant structure. 



V. CONCLUSIONS 

Putative global minima for the Cgo intermolecular po- 
tential of Pacheco and Prates-Ramalho have been lo- 
cated for all clusters containing up to 105 molecules. For 
iV < 15 the structures follow an icosahedral growth se- 
quence, but above this size the global minima are ei- 
ther decahedral or close-packed, with close-packed clus- 
ters becoming more common as the size increases. This 
progression to structures with lower strain energy as the 
size increases is expected, but is more rapid than for most 
atomic clusters, because the PPR potential is much more 
short-ranged than typical interatomic potentials. 

The correspondence between many of the particularly 
low-energy PPR clusters and the high temperature ex- 
perimental magic numbers adds further support to the 
interpretation that these sizes have particularly low free 
energies because they are particularly low in potential 
energy. This interpretation also implies that the icosahe- 
dral magic numbers seen at lower temperature are kinetic 
in origin. The large interfunnel barriers evident from the 
disconnectivity graphs of (C6o)38 and (050)55 add further 



support to this hypothesis by indicating the difficulty of 
major structural transformations in a growing cluster. 
Moreover, preliminarAi-. results from growth simulations 
confirm this scenario O 

However, there are still discrepancies between the high 
temperature magic numbers and the low-energy model 
Geo clusters. Firstly, the magic number at A^=19, which 
persists up to high temperature, is probably due to the 
double icosahedron. This perhaps suggests that the PPR 
potential is still effectively too short-ranged, but an al- 
ternative explanation is the thermal stabilization of this 
structure at the experimental temperature. 

Secondly, a particular preference for structures that 
are based on the Leary tetrahedron is not reproduced. It 
is possible that these structures are only lowest in free 
energy at high temperature, but this seems unlikely for 
two reasons. Firstly, as the entropy is not so strongly 
size-dependent as the energy, magic numbers associated 
with structures that are entropically preferred are likely 
to be less pronounced. However, in the present case the 
magic numbers associated with the structures based upon 
the Leary tetrahedron are the most prominent. Secondly, 
no such transition is seen for (C6o)98- Hence, it is more 
likely that the discrepancy is due to the potential, and 
again suggests that the PPR potential is effectively too 
short-ranged. For the Morse potential the 98-atom Leary 
tetrahedron is lowest in energy for 6.91 < p < 9.45. How- 
ever, the situation is somewhat more complex. Only at 
a few sizes are structures based on the Leary tetrahe- 
dron actually the global minimum of the Morse poten- 
tial. Therefore, the observed preference cannot simply 
be explained by the effective range. It is most likely to 
be related to the orientational degrees of freedom that 
are neglected in a single-site potential, such as the PPR 
model. 

The thermodynamics calculations also indicate further 
discrepancies between experiment and the PPR model. 
We saw in Section [V that icosahedral structures are ther- 



modynamically favoured at high temperature for certain 
sizes, because of their larger vibrational entropy. How- 
ever, except perhaps for 7V=19, there is no experimental 
support for this scenario. 

Although all-atom potentials have hnpn jsniployed for 
modelling clusters of (Ogo) molecules ,E]0'E3 the double 
icosahedron was found not to be lowest in energy for 
iV=19 and the potentials were not applied in the size 
range relevant to structures based on the Leary tetra- 
hedron. Furthermore, these types of all-atom potential 
are unable to reproduce the low-temperature properties 
of bulk Oeo — the molecules have the-, incorrect orienta- 
tion in the low-temperature crystal.EII Although, at the 
temperatures relevant to experimentp-the molecules are 
expected to be able to rotate freely,E3 it could well be 
that the energetic preferences for structures in which the 
molecules can have the preferred orientations will persist 
up to higher temperature. Therefore, to reproduce the 
preference for Leary tetrahedra it is likely that a poten- 
tial that can give the correct orientations in the crystal is 
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required. There are a number of such potentialsjia 
however this orientational preference has sometimes been 
achieved through the introduction of unceahstic electro- 
static properties for the Cgo molecule.Ej StiU, it would 
be interesting to know if these potentials do favour Lcary 
tetrahedra. However, this would be an extremely chal- 
lenging task computationally, both because of the com- 
plexities of the potentials and the huge number of possi- 
ble orientational isomers for such large clusters. 
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TABLE I: Energies and point groups for the putative global minima for (Ceo) at clusters modelled by the PPR potential. 

The energies of these minima when reoptimizcd without the inclusion of the three-body term {E2) arc also given. 
Structural assignments are given for all clusters with N > 13, where i stands for icosahedral, d for decahedral (d+ 
have an over layer on the {111} faces), f for fee and c for close-packed (but not fee). Fee and close-packed structures 
can be unstrained. 
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-10.556666 


48 


C2V 




-50.835900 


-52.786358 


83 


C2V 


c 


-95.639897 


-99.508789 


14 




i 


-11.010918 


-11.384715 


49 


Cs 


d+ 


-51.926724 


-53.919060 


84 


Czv 


d 


-96.907100 


-100.839577 


15 


C2v 


i 


-12.028948 


-12.436957 


50 


Dah 


c 


-53.354529 


-55.400213 


85 


Cm, 


c 


-98.138499 


-102.126340 


16 


C2V 


d 


-13.017348 


-13.429036 


51 


Cs 


c 


-54.449934 


-56.540311 


86 


Csv 


c 


-99.728365 


-103.787568 


17 


C2V 


d 


-14.088380 


-14.538481 


52 


C2v 


c 


-55.797330 


-57.945711 


87 


Cs 


c 


-100.867536 


-104.983413 


18 


Dsh 


d 


-15.163644 


-15.653781 


53 


Cs 


c 


-56.892907 


-59.085964 


88 


Cs 


c 


-102.173549 


-106.335126 


19 


C2V 


d 


-16.218799 


-16.743185 


54 


C2v 


d 


-58.306782 


-60.570092 


89 


Cs 


d+ 


-103.461268 


-107.695215 


20 


C2V 


d 


-17.274543 


-17.832640 


55 


C2v 


d 


-59.408665 


-61.716335 


90 


Cs 


, _|_ 

d+ 


-104.807341 


-109.095407 


21 


Cs 


T _|_ 


-18.353074 


-18.949300 


56 


Cs 


d 


-60.752438 


-63.116143 


91 


Csv 


c 


-106.302975 


-110.678677 


22 


Ci 


T-l- 


-19.425310 


-20.059762 


57 


C2V 


d 


-62.098302 


-64.517510 


92 


Cs 


c 


-107.397152 


-111.806846 


23 




d+ 


-20.686563 


-21.386861 


58 


Dzh 


1 _|_ 

d+ 


-63.386632 


-65.901565 


93 


Cs 


c 


-108.747951 


-113.216329 


24 


Ci 


1 _|_ 

d+ 


-21.759634 


-22.497858 


59 


Td 


c 


-64.769800 


-67.323349 


94 


Cs 


c 


-110.100130 


-114.626428 


25 




1 _|_ 

d+ 


-22.957151 


-23.760958 


60 


Cs 


d 


-65.887129 


-68.480038 


95 


C2V 


d 


-111.653413 


-116.285755 


26 


C2V 


1 _|_ 

d+ 


-24.141029 


-24.968567 


61 


Csv 


f 


-67.169082 


-69.798058 


96 


Cs 


d 


-112.771588 


-117.460346 


27 


Cs 


d+ 


-25.218788 


-26.085977 


62 


Ci 


d 


-68.330220 


-71.021579 


97 


Ci 


d 


-114.123353 


-118.862596 


28 


Cs 


d+ 


-26.403185 


-27.334120 


63 


Cs 


d 


-69.702863 


-72.458970 


98 


Cs 


c 


-115.518307 


-120.295415 


29 


D^h 


d 


-27.542310 


-28.503488 


64 


C2v 


d 


-71.278073 


-74.107991 


99 


Cs 


d 


-117.079750 


-121.959509 


30 


C2v 


d+ 


-28.649061 


-29.682832 


65 


C2v 


d 


-72.382930 


-75.268852 


100 


Td 


c 


-118.491490 


-123.409546 


31 


C211 


d 


-29.963558 


-31.028821 


66 


Cs 


d 


-73.720048 


-76.664850 


101 


Dsh 


d 


-120.052347 


-125.073354 


32 


Ci 


d 


-31.047706 


-32.157893 


67 


C2v 


d 


-75.059516 


-78.062777 


102 


C2V 


d 


-121.170243 


-126.247600 


33 


C2v 


d 


-32.386698 


-33.555532 


68 


Ci 


d 


-76.157881 


-79.200432 


103 


Cs 


d 


-122.516174 


-127.652314 


34 


Ci 


d 


-33.469737 


-34.683456 


69 


Ci 


d 


-77.504125 


-80.605732 


104 


C2v 


d 


-123.864106 


-129.058381 


35 


C2v 


d 


-34.807885 


-36.080228 


70 


Cs 


d 


-78.872607 


-82.039199 


105 


Cs 


d 


-124.981491 


-130.232109 


36 


Ci 


d 


-35.889548 


-37.206754 


71 


C2V 


d 


-80.450302 


-83.690481 












37 


C2V 


d 


-37.226573 


-38.602386 


72 


C7i 


d 


-81.545223 


-84.828030 













